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Abstract 



The stability of five finite difi^erence-time domain (FD-TD) sche- 
mes coupling Maxwell equations to Debye or Lorentz models have been 
analyzed in where numerical evidence for specific media have been 
used. We use von Neumann analysis to give necessary and sufficient 
stability conditions for these schemes for any medium, in accordance 
, with the partial results of 

o\ . 

^ I Keywords : Stability analysis, Maxwell-Debye, Maxwell-Lorentz. 

1/^ '. 1 Introduction 

O 

r"| ■ To describe the propagation of an electromagnetic wave through a dispersive 

medium some extensions to Maxwell equations are used. They involve time 
^ ' differential equations which accounts for the constitutive laws of the material 

that link the displacement D to the electric field E or equivalently the 
polarization P to E. We focus on two of these models (Debye and Lorentz 
^ ' models) which are addressed in [J in view of specific applications to the 

■ interaction of an electromagnetic wave with a human body. In contrast we 

treat any medium which is described by these models. We only consider the 
stability analysis of numerical schemes whereas ^ also treated phase error 
issues. 



*B. Bidegaray-Fesquet is with the LMC-IMAG, CNRS UMR 5523, B.P. 53, 38041 
Grenoble Cedex 9, France. E-mail: brigitte.bidegaray@imag.fr . 
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1.1 Maxwell— Debye and Maxwell— Lorentz models 

In our context (no magnetization) the Maxwell equations read 

(Faraday) 9fB(t,x) = — curlE(t,x), 
(Ampere) atD(t,x) = — curlB(i,x), 

where x G together with a linear constitutive law 

D(t,x) =eoeooE(t,x) +eo / E(t - r,x)x(T)dT, (2) 

where Soo is the relative infinite frequency permittivity and x is the linear 
susceptibility. The discretization of the integral expression Q leads to re- 
cursive schemes (see e.g. [2], jSl)- However, differentiating Eq. ^ leads to 
a time differential equation for D which depends on the specific form of x- 
For a Debye medium 

UdtB + D = trEoeoodtE + eoEsE, (3) 

where tr > is the relaxation time and Eg > £oo is the relative static per- 
mittivity. Defining the polarization by P(t,x) = D(i,x) — eoeooE(t,x), an 
equivalent form is 

t,dtP + P = eo{es-e^)B. (4) 
For a Lorentz medium with one resonant frequency loi, we likewise have 

dfU + udtB + LojU = EQEoodt^ + eoEooi^StE + eoegW^E, (5) 

where > is a damping coefficient, and 

d'fV + udtP + = eo(es - eoo)c^?E. (6) 

If we denote by J the time derivative of P, system (Q) can be cast as 

5tB(i,x) = -curlE(t,x), 

1 (7) 
£o^oo9tE(t,x) = — curl B(t,x) — J(t,x). ^ ' 

^0 
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1.2 Numerical schemes 



A classical and very efficient way to compute the Maxwell equations is the 
Yee scheme We restrict our study to existing Yee based schemes. Other 
methods may be found in the literature in the context of Maxwell-Debye 
and Maxwell-Lorentz equations: see e.g. for pseudo-spectral schemes or 
[H] for finite element-time domain (FE-TD) schemes. 

The Yee scheme consists in discretizing E and B on staggered grids in 
space and time. This allows to use only centered discrete differential op- 
erators. We denote by h the space step (supposed here to be the same in 
all directions in the case of multi-dimensional equations) and by k the time 
step. In space dimension 1, we only consider the dependence in the space 
variable z and classically two polarizations for the field may be decoupled. 
For example, the transverse electric polarization only involves E = and 
B = By. The discretized variables are i?" ~ E{nk,jh) (and similar nota- 
tions for D = Dx) and B^^f ~ B{{n + h)k,{j + and the Yee scheme 
for system (0) reads 

- m) = -(5 7 - 5 7). 



Similarly the Yee scheme for system (jTj) reads 
k ■' ■' noh 3+2 



Usual Maxwell equations consist in taking J- ^ = in Eq. Q or equiv- 
alently 5" = eo^oo-E'" in Eq. (jSJ and leads to a stable second order scheme 
under a Courant-Friedrichs-Lewy (CFL) stability condition. Namely, if 
Coo = denotes the infinite frequency light speed, the CFL con- 

dition reads Cook < h if the space dimension is = 1 and Cook < /i/\/2 for 
A = 2 or 3. 

In contrast to the recursive schemes, we are interested in direct integra- 
tion schemes which are based on the finite difference-time domain (FD-TD) 
discretization of Eqs ((SI) to ® (see 0, 0, 0)- 
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1.3 Outline 



The von Neumann stability analysis is recalled in Sect. |21 We also describe 
the sketch of our proofs which is common for all the schemes. In Section |21 
two one dimensional direct integration schemes for Debye media are pre- 
sented and analyzed, pointing carefully out the physical properties needed 
to ensure stability and the specific cases which have to be handled sepa- 
rately. Numerical applications to physical media are also given. The same 
point of view is carried out for Lorentz media in Sectional Two-dimensional 
results are given in Section [31 

2 Principles of the von Neumann analysis 

The von Neumann analysis allows to localize roots of certain classes of poly- 
nomials, which proves to be crucial here. We recall the main principles of 
this technique. Details and proofs of theorems may be found in pi7] . 

2.1 Schur and von Neumann polynomials 

We define two families of polynomials: Schur polynomials and simple von 
Neumann polynomials. 

Definition 1 A polynomial is a Schur polynomial if all its roots, r, satisfy 
\r\ < 1. 

Definition 2 A polynomial is a simple von Neumann polynomial if all its 
roots, r, lie on the unit disk (\r\ < Ij and its roots on the unit circle are 
simple roots. 

If a polynomial is of high degree or has sophisticated coefficients, it may be 
difficult to locate its roots. However, there is a way to split this difficult 
problem into many simpler ones. For this aim, we construct a sequence of 
polynomials of decreasing degree. Let (p be written as 

(j){z) = Co + ClZ -\ 1- CpZ^, 

where cq, ci . . . , Cp G C and Cp ^ 0. We define its conjugate polynomial cj)* 
by 

4>*{z) = c; + c;_^z + --- + cizP. 

Given a polynomial </>o, we may define a sequence of polynomials 

9m+l[Z) = . 
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It is clear that deg0m+i < d^E^Pm, if (pm ^ 0. Besides, we have the two 
following theorems. 

Theorem 1 A polynomial (pm is a Schur polynomial of exact degree d if 
and only if (prn+i is a Schur polynomial of exact degree d — 1 and |(/>m(0)| < 

IC(o)|. 

Theorem 2 A polynomial (pm is a simple von Neumann polynomial if and 
only if 

• (pm+i is a simple von Neumann polynomial and 1 0^(0)1 < |</>^(0)|, 

or 

• (pm+i is identically zero and 0^ is a Schur polynomial. 

The main ingredient in the proof of both theorems is the Rouche theorem 
(see ^Oj). To analyze 4)q, at each step m, conditions should be checked 
(leading coefficient is non-zero, |(/>m(0)| < |</'m(0)l) ••■) until a definitive 
negative answer arises or the degree is 1. 

2.2 Stability analysis 

The models we deal with are linear models. They may therefore be analyzed 
in the frequency domain. Thus we assume that the scheme handles a variable 
with spatial dependence 

C^" = C/"exp(i^.j), 

where ^ and j € M-^, N = 1,2,3. The amplification matrix G is the matrix 
such that = GC/". We assume that G does not depend on time or on 

h and k separately but only on the ratio h/k. Let (pQ be the characteristic 
polynomial of G, then we have a sufficient stability condition. 

Theorem 3 A sufficient stability condition is that 4>o be a simple von Neu- 
mann polynomial. 

This condition is not necessary. A scheme is stable if and only if the sequence 
(C/")„gN is bounded. Since we assume that G does not depend on time, = 
G^U^ and stability is also the boundedness of (G")„gN- If the eigenvalues of 
G, i.e. the roots r of (/>o, lie inside the unit circle (|r| < 1), then lim^^oo 

G" = 

and the sequence is bounded. If any root lies outside the unit circle then 
G" grows exponentially and the scheme is unstable. The intermediate case 
when some roots may be on the unit circle (and the others inside) may 
lead to different situations. The good case is for example given when G 
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is the identity. Then [/" = U and the scheme is clearly stable. However 
there are other examples of matrices with multiple roots on the unit circle 
that lead either to bounded or unbounded sequences {G'^)nen- We will call 
this property G^-boundedness in the sequel. It is clearly a property of the 
amplification matrix and not of its characteristic polynomial. If the minimal 
stable subspaces associated to the multiple root are one-dimensional then 
is bounded (identity example). If the minimal stable subspaces are 
multidimensional then G'^ grows linearly. Such cases (which occur for our 
schemes) should therefore be handled specifically. 

2.3 Sketch of proofs 

In the next sections, we will not give the proofs, but only list in a table the 
arguments used for each situation. We describe here the general plan and 
give names to specific final arguments used. The detailed proofs may be 
found in for space dimensions 1 and 2. The three dimensional case is 
much more tedious and is work in progress. 

Usually the system is given in a implicit form. The first step consists 
in writing it in an explicit form. This yields the amplification matrix G. 
Then we compute its characteristic polynomial (f)Q. In order to perform a 
von Neumann analysis, we compute the series {4>m)- In the general case, 
under the assumption that the stability condition cannot be better than 
Maxwell's, we can apply either Theorem ^ ( J'/ieorem Q argument) or The- 
orem |^ ( T/ieorem argument), check estimates at each level until (j)m is a 
one degree polynomial. Special cases arise when Eg = £oo, sin(^/2) = or 
±1, and sometimes for limit values of physical coefficients. In these cases, 
different points of view have to be considered: 

• Theorem 121 has to be used instead of Theorem^ 

• Some eigenvalues lie on the unit circle (mostly ±1 or zizi) and are 
simple, it is then sufficient to study only the other eigenvalues {sub- 
polynomial argument) and we conclude to a simple von Neumann poly- 
nomial and stability, 

• Some eigenvalues lie on the unit circle and are not simple, and besides 
the study of the other eigenvalues (to prove that the polynomial is 
a von Neumann one), we have to find out if the associated minimal 
stable subspaces are one- (stable case) or multidimensional (unstable 
case). This may be checked directly on the form of matrix G {G form 
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argument), or necessitates the computation of eigenvectors {eigenvec- 
tors argument). If only one eigendirection is found for a multiple 
eigenvalue, the minimal subspace is necessarily multidimensional. 



3 Debye media 

We address two discretizations of Maxwell-Debye equations. The first one 
uses a (B,E,D) setting for the equations and the second a (B,E,P,J) 
formulation. 

3.1 Debye-Joseph et al. model 

In [S], Joseph et al. close System © by a discretization for Eq. @, namely 



eoEooir ; 1- eoEs — ir ; 1- 



k " ° 2 k 2 

(10) 

System (|5|)- (fTm) may be cast in an explicit form which handles the variable 

t/; = \cooB'^-lElD]/eoeoo) 
and the amplification matrix G reads 

/ 1 -A(e*«-1) 

(l+^)A(l-e-'?) (l-feO+(l+^)A^(e'g-2+e-'g) 2<5 

\ -A(l-e-^«) A2(e^« -2 + e-^«) 1 

where A = Cook/h is the CFL constant, 6 = k/2tr > is the normalized 
time step and e'^ = Es/Soo > 1 denotes the normalized static permittivity. 
Moreover we define 

q = -X\e'^ - 2 + e-'^) = AX'^ sm\^/2). 

The characteristic polynomial is proportional to 

MZ) = [l + Sei)]Z^-[3 + 6ei-il + 5)q]Z^ 
+ [3-6ei-{l-6)q]Z-[l-6ei]. 

The proofs are summed up in Table ^ and we deduce that the stability 
condition is q < 4 if Ss > £oo and g < 4 if eg = Soo • 
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q 




argument 


result 


]0,4[ 


> Eoo 


Theorem ^ 


stable 


]0,4[ 


— £oo 


Theorem |21 


stable 





> eoo 


G form 


stable 


4 


> eoo 


Theorem |2l 


stable 


4 




eigenvectors 


unstable 



Table 1: Proof arguments and results for the Debye- Joseph et al. model. 
3.2 Debye- Young model 

In [U], Young closes System ® by two discretizations for Eq. Q, namely 
tr^ r^— = + eo{e, - Eoo)E^, (11) 



= + eo(e. - eoo)^^ (12) 

n+i 

Although J- ^ is used for the computations, this not a genuine variable for 
System (|9])- (|TT|) ~ (fT2|) which handles the variable 

1 1 

rn t 



and the amplification matrix G reads 

/ 1 -A(e^€-1) 

A(l-e-'S) l+5~5a+35'^a-{l+S)q is 25 

l+5a (l+5)(l+(5«) l+(5 l+(5« 

n 25a 1^ 

\ ^ 1+5 1+5 

with the same notation as above and a = — 1 > 0. 
The characteristic polynomial is proportional to 

MZ) = [{l + Sa){l+6)]Z^ -[3 + 6 + 6a + 36'^a-{l + 5)q\Z^ 
+ [3-6-6a + - (1 - S)q]Z - [(1 - 6a){l - 6)]. 

Again, the proofs are summed up in Tabled 

The stability condition is therefore g < 4 and 5 < 1 if Eg > £oo and q < 4 

if Es = Eoo- 
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q 


£s 


6 


argument 


result 


]0,4] 


> £00 


]o,i[ 


Theorem ^ 


stable 


]0,4[ 


— £^00 


> 


Theorem 12 


stable 





> £00 


> 


G form 


stable 


]0,4] 


> £00 


1 


sub-polynomial 


stable 


4 


£00 


> 


eigenvectors 


unstable 



Table 2: Proof arguments and results for the Debye- Young model. 

3.3 Conclusion for one-dimensional Debye schemes 

If Eg > Soo, the pure CFL condition g < 4 is the same for both models. 
It is exactly the condition for Maxwell equations. However Young model 
necessitates another condition, 5 < 1, which corresponds to a sufficient dis- 
cretization of Debye equation (jlj . Even if we are interested here in stability 
properties, such conditions are to be taken to ensure equations to be cor- 
rectly taken into account. Results are given in physical variables in Table 



Scheme 




dimension 1 




Joseph et al. 


q<4 




Young 


g < 4, 6<1 


k < min(-^,2tr) 

^ Coo 




Joseph et al. 


q < 4 




Young 


g < 4 


k<^ 



Table 3: Stability of Debye models for Es > £00 and £s = £00 ■ 

To compare conditions on q and 5, let us consider a simple physical 
case. We assume that a matter with Soo = 1 (and thus Cqo — SlO^ms"-*^) 
is lighted by an optical wave of say wavelength 1 ^um. The space step h 
has to be smaller than this wavelength, and therefore q < 4 reads at least 
A; < 10~^^ s. In a Debye medium, relaxation times are of the order 
of a picosecond (or even a nanosecond) which is many decades larger than 
the previous bound. The estimate q < 4 is thus predominant and both 
models present the same advantages. Only the value of £00 yields the CFL 
condition. A typical example is water for which = 

1.8, es = 81.0 and 

tr = 9.4 10"^2 s 0. Condition k < 2t^ comes to A; < 1.88 10"^^ s. Condition 
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g < 4 yields a similar condition if /i = 4.2 10~^ m. This is of course much 
larger than any reasonable space step for Maxwell equations and optical 
waves. The stability condition for water is q < 4 for both schemes. A quite 
different material is for example the 0.25-dB loaded foam given in for 
which Eoo = 1.01, ffs = 1-16 and tj. = 6.497 10-^° s. Condition k < 2tj. comes 
to k < 1.310~^s and g < 4 yields a similar condition if h = 3.910~^m. 
Once more, the stability condition for water is q < 4 for both schemes. 

In conclusion for current material the stability condition is the same 
for Maxwell-Debye equations as for the usual Yee scheme. The result an- 
nounced in 1 was q < 4 for Joseph et al. scheme and for water, which is 
consistent with our result. 



4 Lorentz media 

Three discretizations of Maxwell-Lorentz equations are now addressed. The 
first one uses a (B, E, D) setting and the two others a (B, E, P, J) formula- 
tion, but differ from the time-discretization of J. 

Each of these models reads the same in the harmonic (z^ = 0) or an- 
harmonic {v > 0) cases. However the analysis will differ greatly since (pi = 
for all the schemes in the harmonic cases. 



4.1 Lorentz— Joseph et al. model 

In [S], system Q is closed by a discretization for Eq. Q, namely 



^„+l 


- 2E^ + EJ-^ 




^2 

- 21?; + U^-^ 


A;2 



VEqEc 



+ ly 



2k 

'7 



2k 



2^] 



n+1 



+ D 



n-1 ■ 



(13) 

The explicit version of system (|8|)- (|13j) does not use explicitly the value of 
DJ^^ and therefore this system handles the variable 



\c^B'^J^,E^,E^-\DJ/eoea 



The amplification matrix G reads 



V -A(l 



1 

2.5A(l-e-'^) 





-A(e^« - 1) 

2-q{l+5+uj) 
1+5+ujei, 
1 












2uj 



1 
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where 6 = vk/2 > is the new normahzed time step, and uj = lo\1? 12 > 
denotes the normahzed squared frequency. The other notations used for the 
Debye model remain vahd. 

The characteristic polynomial is proportional to 

(/>o(Z) = [l + 5 + a;4]Z^-[4 + 2(^ + 2t^4-(l + 5 + w)g]z3 

+ [6 + 2uje!^ - 2q]Z'^ -[A -25 + 2Loe'^ -{1-5 + uj)q]Z 
+ [l-5 + u;e'^. 

The proofs are summed up in Table |1] for the an-harmonic and the harmonic 
case. 



Q 




argument 


result 


an-harmonic: u > 


]0,2[ 


> £oo 


Theorem ^ 


stable 


]0,2] 


— ^oo 


Theorem 21 


stable 





> £00 


G form 


stable 


2 


> £00 


sub-polynomial 


stable 


harmonic: = 


]0,2[ 


> £00 


Theorem |21 


stable 


]0,2] 




sub-polynomial 


unstable 





> £00 


G form 


stable 


2 


> £00 


sub-polynomial 


stable 



Table 4: Proof arguments and results for the Lorentz- Joseph et al. model. 

In the an-harmonic case the stability condition is g < 2 whatever Eg > 
is. The Es = £00 harmonic case, needs some explanation. For q s]0, 2], 
may be cast as the product of two second order polynomials. The roots 
are two couples of conjugate complex roots of modulus 1. For the specific 
value q = 2uj/{\ + uj), which always lies in the interval ]0, 2], the two couples 
degenerate in one double couple, and the associated minimal stable sub- 
spaces are two-dimensional. To avoid this instability one may think to bound 
q and say that the scheme is stable provided q £ [0, 2a;/(l + uj)[. But if we 
come back to the original variables, we see that this is not an upper bound 
on k but rather a lower bound on /i, which we surely do not want. It is 
therefore better to avoid using Joseph et al. scheme in this very specific 
case, = £00 and = 0, and we hope to find a better scheme for this case 
in the following examples. 
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4.2 Lorentz-Kashiwa et al. model 

In |7], Kashiwa et al. close a modified version of System Q, which consists 
of the three first equations in System ((TUl . by a discretization for Eq. 
namely 



£0£oo /^n+l 

k ^ 
1 



1 



1 



1 



_(pn+l 
k^ i 



The explicit version of system p4|) handles the variable 



(14) 



and the amplification matrix G reads 

/ 1 -A(e*« - 1) 

-A(l-e-*«){A-|a;«) A-gA-(2-g)|(^a 



V 



A 



A 

{2—q)LL>a 
A 





LP 

A 



\ 



-1 

A 

\-w l_ 
A A 

-2a; 2-A 
A A 



where together with the previously defined notations, A = 1 + 5 + oje'g/2. 
The characteristic polynomial is proportional to 

MZ) = [l + <5 + la;4]Z^-[4 + 25-(l + 5 + ia.)(?]z3 

+ [6 - uje', + - 2)g]z2 - [4 - 25 - (1 - 5 + \u:)q]Z 



+ [l-5+-^e'J. 

The proofs are summed up in Table |SJ Both in the an-harmonic and har- 
monic cases, the stability condition is g < 4 which is much better than the 
previous scheme since we gain a factor 2 on A; and we have no problem when 
Es = Eoo and = as for the previous model. 
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q 




argument 


result 


an-harmonic: i/ > 


]0,4[ 


> Boo 


Theorem ^ 


stable 


]0,4[ 


— ^oo 


Theorem |21 


stable 





> £00 


G form 


stable 


4 


> £00 


eigenvectors 


unstable 


harmonic: u = Q 


]0,4[ 


> £00 


Theorem |21 


stable 





> £00 


G form 


stable 


4 


> £00 


eigenvectors 


unstable 



Table 5: Proof arguments and results for the Lorentz-Kashiwa et al. model. 



4.3 Lorentz— Young model 

In [S], System Q is closed by a discretization for Eq. ©, namely 



n+l 



p") = r 



2 



(15) 



+u;?(es-£oo)£o^"-c^?^". 
The explicit version of System (|9|)- (|15|) handles once more the variable 

[/; = *(cooi?j;| , E], P;/eo£oo, A:j;/£0£oo) 
and the amplification matrix G reads 



-A(e*« - 


1) 








\ 


(l-g)(l+<5)- 






1-5 




1+5 




1+5 


1+5 








l+5-2i^ 


1-5 




l+<5 




1+5 


1+5 










1-5 


/ 


l+<5 




1+5 


1+5 



/ 1 

-A(l-e- 


V 

The characteristic polynomial is proportional to 

(^o{Z) = [l + 5]Z^ -[^ + 25 -2uje',- {I + 5)q]Z^ 
+2[3 - 2Loe'^ + (cj - l)q]Z'^ 
-[4 -25- 2uje'^ - (1 - 5)q]Z + [1 - 6]. 
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The proofs are summed up in Table IHl This scheme combines three draw- 
backs we have aheady encountered. First as for the Debye model, there 
is an extra condition on the time step: oj < 2/(2eg — 1). This will have 
to be compared to the condition on q for physical examples. Second, as 
for the Lorentz-Joseph et al. scheme we need a twice smaller k than for 
raw Maxwell equations: q < 2 instead of g < 4. Last, and also as for the 
Lorentz-Joseph et al. model, the Es = £oo and v = {) leads to an instability. 
This is exactly the same story. This time q = 2lo leads to double couples of 
conjugate complex roots of modulus 1, with two-dimensional minimal stable 
sub-spaces. If w < 1 this value of q is however never reached, but a; < 1 is a 
stronger assumption than to < 2/{2£'^ — 1). We will see what this amounts 
to in numerical applications. 



q 


Es 


UJ 


argument 


result 


an-harmonic: u > 


]0,2[ 


> Eoo 


< ^ 


Theorem ^ 


stable 


2 


> Eoo 




]0,2] 


— ^oo 


< 2 


Theorem |21 


stable 


]0,2] 


— '^oo 


= 2 


sub-polynomial 


stable 


2 


> £oo 


'2 


Theorem |2 


stable 





> £oo 


< ^ 


G form 


stable 


harmonic: = 


]0,2[ 


> £oo 


< ^ 


Theorem |21 


stable 


2 


> Eoo 




]0,2] 


— ^oo 


< 2 


eigenvectors 


unstable 


]0,2] 


— ^oo 


= 2 


Theorem [2 


stable 


2 


> £oo 


2 

2e',-l 


eigenvectors 


unstable 





> £oo 


- 2e^-l 


G form 


stable 





— ^oo 


< ^ 

^ 2e^-l 





— ^oo 


2 

~ 2e^-l 


eigenvectors 


unstable 



Table 6: Proof arguments and results for the Lorentz- Young model. 

4.4 Conclusion for one-dimensional Lorentz schemes 

We can summarize all our results for Lorentz schemes in Table [7| We chose 
not to translate the result for the Young scheme for condition 
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on h (q < 2uj) but as a condition on k (co < 1, and therefore g = 2a; is not 
reached) . 

For the harmonic Young scheme if Ss > £00 the condition is slightly better 
since q = 2 and uo < 2/(24 ~ 1)' or g < 2 and uj = 2/(24 ~ 1) ^1^° yield 
stable schemes. 

Contrarily to Debye materials, for which Joseph et al. model and Young 
model compete, the Kashiwa et al. model seems to overcome others for 
Lorentz material. First, there is a gain in CFL condition q < 4 is twice 
better as q < 2, second, there are no instabilities for limiting values of 
the physical coefficients and last there are no extra condition on the time 
step. In practice, an extra condition is however needed to account for the 
dynamics of the Lorentz equation, but not for stability reasons. 

However we can compare the relative strength of the different conditions 
on k for Joseph et al. and Young models. The values used in jjj are £00 = 1, 
Es = 2.25, ui = 410i*5rads"^ and ly = 0.56 10^*^ rads"^ Condition lo < 
2/y/24 - 1 comes to A; < 2.7 lO"^"^ s which is very small and corresponds to 
h = 1.13 1"^ m in the q < 2 condition. This space step is more than sufficient 
to discretize optical waves. For such a material the extra condition imposed 
by the Joseph et al. scheme is stronger than the basic CFL condition. The 
Kashiwa et al. model is then more advisable. 

In jni there is a totally different material for which £^0 = 1.5, Eg = 3, wi = 
27r 5 10^*^ rads~^ and u = lO^'^rads"^ (these round values certainly refer to 
a model material). In this case io < 2/ 2eg — 1 comes to k < 3.610~^^s 
which corresponds to h = 1.9 1~^ m in the q < 2 condition. For this material 
condition q < 2 is the strongest for optical waves. The Kashiwa et al. model 
is however more advisable, since it allows g < 4 instead of g < 2. 

The results obtained in [T] where obtained for our first cited material 
and for Joseph et al. and Kashiwa et al. models. He observed instabilities 
for ^ > f . We note that if C < f then sin(^/2) < 1/^2 and g < 2 instead of 
g < 4. This is exactly our result. He found also the Kashiwa et al. scheme 
to stable for g < 4. 

5 Two-dimensional results 

In a two-dimensional context where unknowns depend only on space vari- 
ables X and y, Maxwell system may be split in two decoupled systems cor- 
responding to the transverse electric (TE) {B^, By, E^) and the transverse 
magnetic (TM) {Bz, Ex, Ey) polarizations. In the one-dimensional case. 
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Scheme 




dimension 1 


an-harmonic: v > 0, and Eg > £00 


Joseph 


q<2 




Kashiwa 


g < 4 


Coo 


Young 


<?<2, 


k < min(^|— , 


harmonic: 1^ = 0, and Eg > £00 


Joseph 


q<2 




Kashiwa 


g < 4 


(-00 


Young 


q<2, 
^ < 2el-l 


< min(^|— , — ) 


harmonic: 1^ = 0, and £3 = £00 


Joseph 


a < — 


condition on /t 


Kashiwa 


g < 4 


k<^ 

(--00 


Young 


q<2, 

LO <1 





Table 7: Stabihty of an-harmonic and harmonic Lorentz models for £s > £00 
and £s > £oo- 

Maxwell-Debye equations were represented by three equations and Maxwell- 
Lorentz by four equations. In the TE polarization, one more Faraday equa- 
tion is added and we have four equations for Maxwell-Debye and five equa- 
tions for Maxwcll-Lorcntz. In the TM polarization for the Maxwell-Debye 
model, one Ampere equation and one Debye equation have to be added, 
leading to five equations systems. For the Maxwell-Lorentz model, there 
are one Ampere equation and two Lorentz equations more, and the system 
consists of seven equations. 

The principle of the stability analysis is exactly the same, but we now 
have larger polynomials to study. A small miracle however happens: one- 
dimensional polynomials are a factor in two-dimensional polynomials. More 
precisely we now denote by and hy the space steps in the x- and y- 
directions respectively and by q the quantity 

/ ^2 \ 

q = qx+qy = 4cl^(-r2 sm\^^/2) + sm^i^y/2) 
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(recall q = 4c^^ sm'^{^x/'^) in ID)- Then in the two-dimensional TE polar- 
ization 

0f'^^(Z) = [Z-l],^i^(Z), 

for all the Maxwell-Debye and Maxwell-Lorentz schemes we study here. 
This could be a problem, if 1 is already a root of (^g^(Z), i.e. when q = 0, 
but it happens that it is never a problem: minimal stable sub-spaces are 
always one-dimensional. In the TM polarization, the same factorization 
occurs but the remaining polynomial is slightly more complicated, namely 

c/>f'™(Z) = [Z-l]^o(^)0^^(^), 

where ipo{Z) is equal to: 

- Debye-Joseph et al. model 

[(l + &^)Z-(l-fe'J]. 

- Debye- Young model 

[(1 + q)(1 + 6a)Z - (1 - a)(l - 5a)]. 

- Lorentz-Joseph et al. model 

[{1 + 5 + uje',)Z^ -2Z + {l-5 + toe',)]. 

- Lorentz-Kashiwa et al. model 

[{1 + 5 + \i^e[)Z^ - (2 - ^4)Z + {l-5 + 

- Lorentz- Young model 

[(1 + <5)Z2 - 2(1 - we'jZ + (1 - 5)]. 

As for the TE polarization the extra eigenvalue 1 is never a source of in- 
stability. The other extra eigenvalues always lie inside or on the unit circle 
(conjugate complex roots). The only problem is when modulus 1 eigenvalues 
are also eigenvalues of the one-dimensional polynomial. This only occurs for 
the Lorentz-Joseph et al. scheme is Es = £00, and q = 2a;/ (1 +0;), which is 
a resonant value we have already pointed out in the harmonic case for this 
scheme. 

We shall not duplicate Tables 01 and [3 for two-dimensional models. If 
hx = hy = h, condition g < 4 becomes k < /i/(\/2coo) and condition q <2 
becomes k < h/{2coo) in the physical variables. Besides, Lorentz-Joseph et 
al. model which was leading to a lower bound on h in the harmonic case, 
leads also to such a bound in the an-harmonic case. These are the only 
differences with Tables IHl and |3 
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6 Conclusion 



We have studied a class of FD-TD schemes for dispersive materials based on 
the Yee scheme for Maxwell equations and compared them from the stabil- 
ity point of view. This study was inspired by Petropoulos ; 1_ who performs 
the same analysis but using specific values for the physical and numerical 
constants and using numeric routines to locate eigenvalues of the amplifi- 
cation matrix. Here we have general results which gives you the constraint 
on numerical constants {k and h) for any Debye or Lorentz material. Our 
results confirm those of Petropoulos. 

For usual Debye media, both studied schemes are stable under the same 
conditions as the Yee scheme, ensuring also, if applied to optical waves, a 
fine discretization of the Debye equation. Among the studied schemes for 
Lorentz media, Kashiwa et al. model clearly ranks first as far as stability is 
concerned., Its stability condition is also that of the Yee scheme. However 
to take properly into account the Lorentz model, a smaller time step may 
have to be chosen, independently of stability issues. Such results have been 
proved for ID and 2D models. The 3D case, which is much more tedious, is 
being studied and analogous results are expected. 
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